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Random constraint satisfaction problems are interesting model systems for spin-glasses and glassy 
dynamics studies. As the constraint density of such a system reaches certain threshold value, its 
solution space may split into extremely many clusters. In this work we argue that this ergodicity- 
breaking transition is preceded by a homogeneity-breaking transition. For random A-SAT and 
A'-XORSAT, we show that many solution communities start to form in the solution space as the 
constraint density reaches a critical value a cm , with each community containing a set of solutions 
that are more similar with each other than with the outsider solutions. At a cm the solution space is 
in a critical state. The connection of these results to the onset of dynamical heterogeneity in lattice 
glass models is discussed. 

I. INTRODUCTION 

Random constraint satisfaction problems (CSPs) have attracted a lot of research interest from the statistical physics 
community in recent years^^. These are model systems for understanding typical-case computational complexity of 
nonpolynomial-complcte problems in computer science, some of them have also important applications in modern 
coding systems, such as the low-density-parity-check codes . T he energy landscape of a no ntrivial random CSP problem 
is usually very complicated, similar to those of spin-glasses^' and lattice glass models^'. Therefore understanding the 
configuration space property of random CSPs is also very helpful for developing new insights for spin-glasses, glassy 
dynamics, and the jamming phenomena of colloids and granular systems. 

As the density a of constraints increases, the solution space of a CSP will experience a series of phase transitions^. 
One of them is the clustering transition at certain threshold constraint density ay> where the solution space splits 
into exponentially many isolated solution cluste rs or Gibbs states. This ergodicity-breaking transition has very 
significant consequences for glassy dynamic d 10 * 11 * and stochastic local search processea^Ell. Numerical simulations^! 
further suggested that, before ergodicity of the solution space is broken, the solution space has already been non- 
homogeneous, with the formation of many solution communities. In this paper we determine the critical constraint 
density a cm for the solution space of a random CSP to become heterogeneous. We find that a cm < ay, and that at 
d = oi crn the solution space is in a critical state, in which the boundaries between different solution communities of 
the solution space disappears, while at a > a cm the solution space contains many well-formed solution communities. 

Heterogeneity of the configuration space of a complex system can cause heterogeneity in the dynamics of this system. 
The results of this work may be helpful for understanding more quantitatively the nature of dynamical heterogeneity 
in supercooled liquids^HUl an d lattice glass models^'. 

II. THEORY 

A constraint satisfaction formula has N vertices k, . . .) and M constraints (a, 6, c, . . .), with constraint density 
a = M/N. A configuration of the model is denoted by a = {ci, 02, . . . , ct/v}, where <Tj = ±1 is the spin state of vertex 
i. Each constraint a represents a multi-spin interaction among a subset (denoted as da) of vertices, and its energy 
E a is either zero (constraint being satisfied) or positive (unsatisfied). For example, E a may be expressed as 

E a = 1 - J a J] f< (!) 

where J a = +1 or —1 depending on the constraint a. The total energy of a spin configuration a is the sum of 
individual constraint energies, E(a) — 53 =l ^a- 

A solution of a constraint satisfaction formula is a spin configuration of zero total energy. The whole set of solutions 
for a given energy function E(a) is denoted as S and referred to as the solution space. The energy landscape of a 
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FIG. 1: Evolution of the solution space S of a constrained spin system: At low constraint density a (left panel), 5 is 
homogeneous and the solution-pair entropy density s(q) is a concave function of the overlap q. Solution communities start to 
form as a exceeds a threshold value a cm (middle panel); S then becomes heterogeneous and the function s(q) changes to be 
non-concave. An ergodicity- breaking transition occurs as a reaches a larger threshold value «d (right panel), where the solution 
communities separate into different solution clusters and s(q) becomes non-monotonic. As a further increases, solution-pairs 
with intermediate overlap values may disappear completely, and then s(q) is not defined for these intermediate overlap values. 



CSP may be very complex and its ground-state degeneracy usually is extremely high. To investigate the structure of 
a solution space S, we define a partition function as 

JV 

E E exp^E'i ?) , (2) 

s 1 ess 2 es i=i 

where a; is a binding field. Each solution pair (a 1 , a 2 ) contributes a term exp[Nxq(a 1 , a 2 )} to Z(x), with q^a 1 ,^ 2 ) = 
(1/N) Yli=i °i °f being the solution- solution overlap. We introduce an entropy density s(q) by expressing the total 
number Af(q) of solution- pairs with overlap q as Af(q) = exp[Ns(q)]. Then Eq. ^ is re- written as 

Z(z)= ^ exp[iV (*(«?) + . (3) 

The free entropy^ $(cc) of the system is related to the partition function by = \nZ(x) In the thermodynamic 
limit of N —> oo, the free entropy density 4>(x) = $>(x)/N is related to the entropy density s(q) by 

<p(x) = max s(q) + xq = s(q(x)) + xq(x) . (4) 
«e[-i,i]L -I 

In Eq. Q, q(x) is the overlap value at which the function s(q) + xq achieves the global maximal value. q(x) is also 
the mean value of solution-solution overlaps under the binding field x. 

At x = 0, the maximum of Eq. Q is achieved at q(0) = qo, the most probable solution-pair overlap value; at 
the other limit of x — > oo, q(po) = 1. If the entropy density s(q) is a concave function of q G [qo, 1] (Fig. [I] left 
panel), then for each x > there is only one mean overlap value q, and q(x) changes smoothly with x. On the 
other hand, if s(q) is non-concave in q G [qq,1] (Fig. [TJ middle and right panel), then at certain value x* of the 
binding field, there are two different mean overlap values, and the value of q(x) changes dis continuously at x = x* 
(a field-induced first-order phase-transition). In this work, we exploit this correspondence between the non-concavity 
of s(q) and the discontinuity of q(x) to determine the threshold constraint density a cm at which the solution space 
S becomes heterogeneous. Many solution communities can be identified in a heterogeneous solution space ^ Each 
solution community contains a set of solutions which are more similar with each other than with the solutions of other 
communities. These differences of intra- and inter-community overlap values and the relative sparseness of solutions 
at the boundaries between solution communities cause the non-concavity of s(q). 



III. APPLICATION TO THE RANDOM K-SAT PROBLEM 



We begin with the random K-SAT, a prototypical CSF * 3 * 4 * 9 ( In a random i\-SAT formula, the number of vertices 
in the set da of each constraint a is fixed to K, and these K different vertices are randomly chosen from the whole 
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set of ./V vertices. Depending on the spins of these K vertices, the energy of a constraint a is either zero or unity: 

Ea=W 1 ^, (5) 

where J l a — ±1 with equal probability. The solution space S of a large random if-SAT formula is non-empty if the 
constraint density a is less than a satisfiability threshold a s (Kj^. Before a s (K) is reached, S has an ergodicity- 
breaking transition at a clustering transition point otd(K), where it breaks into extremely many solution clusters^. 
We will see shortly that this clustering transition is preceded by another transition at a = a cm (K) , where S starts to 
be heterogeneous as many solution communities are formed. 

We use the replica-symmetric cavity method of statistical mechanics^ to calculate the mean overlap value q(x) at 
a < o>d{K). As the partition function Eq. Q is a summation over pairs of solutions (a 1 , c? 2 ), the state of each vertex 
is a pair of spins [a, a'). Consider a vertex i which is involved in a constraint a, i € da. The following two cavity 
probabilities pi^ a {&i, and p a ^i{(Ji 1 a[) arc defined: Pi^ a {&i, &[) is the probability that, in the absence of constraint 
a, vertex i has spin value Oi in solution a 1 and value a[ in solution <r 2 ; and p a ^i(ai, a[) is the probability that the 
constraint a is satisfied conditional to vertex i being in state (er^er^)- One can write down the following iterative 
equations: 



j£da\i cr j£da\i & 

+fcf a 5-f a II P^a(-Jl-Ji) , (6) 
j£da\i 

Pi^aicruv'i) = Ce^< J] pb^tfaoi) , (7) 

b£di\a 

where <5^ is the Kronecker symbol, C is a normalization constant, and di denotes the set of constraints that vertex i 
is associated with. The probability Pi((Ji, &[) of vertex being in the spin-pair state (<7i, &[) has the same expression as 
Eq. Q but with di\a replaced by di. In writing down the above cavity equations, we have applied the Bethe-Peierls 
factorization approximation of cavity probabilities, which corresponds to the replica-symmetric cavity theory^M For 
each vertex i the probabilities Pi and Pi^. a have the symmetry that Pi{+, —) = pi(—, +) and p^ a (+, — ) = p^ a (— , +). 
The mean overlap is expressed as 



1 N 

jvX^ (+'+)+«(--) ~ 2 Pi (+'-)] > ( 8 ) 



and the free entropy density <fi(x) can also be expressed by the cavity probabilities^. The overlap susceptibility 
X = dq(x)/dx is a measures of the overlap fluctuations, 



N N 



= 1 0=1 

where (...) means averaging over solution-pairs under the binding field x. 

Equations ^ and Q can be solved by population dynamics methocP. Some of the analytical results as obtained 
for the random 3-SAT problem are shown in Fig. [5] (the results for K > 4 are qualitatively the same). When 
q < a cm = 3.75, the mean overlap q increases with the binding field x smoothly, indicating that the solution space 
S of the random 3-SAT problem is homogeneous. The overlap susceptibility x{ x ) nas a single peak, whose value is 
inverse proportional to (a cm — a) and diverges at a = a cm and x — x cm = 0.0024. The susceptibility x(x) is again 
finite when a exceeds a cm , but the mean overlap q(x) changes discontinuously with x at certain threshold value x*. 
This first-order phase transition at a > a cm suggests that in the space S many solution communities (groups of 
similar solutions) are formed. For x > x* the partition function is predominantly contributed by intra-community 
solution-pairs (overlap favored), while for x < x* it is contributed mainly by inter-community solution-pairs (entropy 
favored). The different solution communities of S all belong to the same solution cluster (s(q) is non-negative for 
any q £ [<7o>l]) as long as a is less than ad — 3.87 9 , but at a = ad they start to break up into different solution 
clusters (s(q) is not defined for some intermediate q values^) . At a = a cm the solution space S is in a critical state 
at which the boundaries between different solution communities disappear. This situation is qualitatively the same 
as the critical state of water at 647K and 22.064MPa, where the liquid and the gas phase are indistinguishable. 
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FIG. 2: The mean overlap q(x) (A) and the overlap susceptibility x{ x ) (B) at different constraint density values a for the 
random 3-SAT problem. In (A) the value of a increases from 3.72 to 3.76 with step size 0.01. The insets of (B) show that the 
peak value of \ diverges inverse linearly with a and x as the critical point (a cm = 3.75, x cm = 0.0024) is approached (dashed 
lines are linear fittings). 




FIG. 3: Mean overlap value of solutions with a reference solution 3* for a single random 3-SAT formula oi N — 10 vertices 
and constraint density a = 3.85. Solid lines are mean-field analytical result on this single formula^, and the symbols with 
error bars are single-spin flips simulation results. Each sampled solution trajectory starts from a* and is equilibrated for at 
least 10 7 Monte Carlo steps (each step corresponds to TV spin-flip attempts). More than 1000 overlap values with a* are then 
sampled at time interval of 10 4 Monte Carlo steps. The blue dashed line marks the equilibrium transition value of x. 



For the random 4-SAT problem, we find that a cm = 8.4746, which is consistent with the simulation results of 
RefPl. The value of OL cm is much below the clustering transition point — 9.38-4 

The solution space heterogeneity can also be detected using single solutions as reference points^. Figure [3] shows 
the theoretical and simulation results on a random 3-SAT formula with N = 10 5 vertices and M = 3.857V constraints. 
The reference solution a* is uniformly randomly sampled from the solution space, and single-spin flips are used in 
the simulation to sample solutions a with weight proportional to exp[Nxq(a* , <j)}^. The replica-symmetric cavity 
method predicts an equilibrium discontinuous change of the mean overlap value with solution a* at x — x* — 0.00267, 
which is confirmed by simulations. For x > x* most of the sampled solutions are in the same solution community 
of <?*, but for x < x* the sampled solutions are scattered in many different solution communities. Because of the 
high degree of structural heterogeneity, at x < x* it takes about 10 7 Monte Carlo steps to travel from one solution 
community to another different community, making it very difficult to sample independent solutions. 

When the solution space of the random iiT-SAT problem becomes heterogeneous at a > a cm (K), the replica- 
symmetric cavity theory, which leads to Eqs. ([6|-([7]), probably is not sufficient to describe its statistical properties. 
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FIG. 4: Schematic phase diagram for a constraint satisfaction problem, using temperature T and constraint density a as control 
parameters. The configuration space is homogeneous and ergodic in region I. As the temperature T decreases to T cm (a), a 
homogeneity-breaking transition occurs, and the configuration space becomes non-homogeneous but still ergodic (region II). 
As T further decreases to Td(a), an ergodicity-breaking (clustering) transition occurs, and the configuration space breaks into 
many separated clusters (region III). At T = 0, the ground-state configuration space is non-homogeneous at a > a cm and 
non-ergodic at a > ad- 

In a future publication we will report the result of the stability analysis on the replica-symmetric cavity equations, 
and present a mean-field study using the first-step replica-symmetry-breaking cavity theory. 

IV. APPLICATION TO THE RANDOM A-XORSAT PROBLEM 

The A'-XORSAT problem has wide-spread applications in low-density-parity-check codeiP and is also extremely 
studiecPSHUl The constraint energy E a of this model is expressed in Eq. Q, where J a = ±1 with equal probability. 
The solution space of a random AT-XORSAT problem breaks into exponential solution clusters of equal size at a 
clustering transition point a = ad(K) 22 ' 23 '. We have applied the replica-symmetric cavity method to this problem 
and obtained the same qualitative results as for the random AT-SAT problem, namely that before the ergodicity of 
the solution is broken, exponentially many solution communities start to form in the solution space as the constraint 
density reaches a critical value a cm (K). At a = a cm (K) the solution space is in a critical state. For K = 3 we find 
that a cm (3) = 0.6182, which is much lower than the value of ad(3) = 0.818^-R For the random 4-XORSAT problem, 
we find that a cm (4) = 0.504, while ay (4) = 0.7722H. 

The random AT-XORSAT problem has a gauge symmetry that can be exploited to simply the mean-field 
calculations^. Suppose a 1 is a solution, we can perform a gauge transformation ai — > crj = Ui<j} to change the 
constraint energy Eq. Q into E a = 1 — Iljeaa^- All the coupling constants J a then become unity. The solution 
space structure of the random AT-SAT problem looks the same from any a reference solution. We have used this nice 
property to calculate the total number of solutions that have a overlap value q with a randomly chosen reference 
solution. 



V. DISCUSSION 

The main conclusion of this work is that, the solution space of a random constraint satisfaction problem has a 
transition to structural heterogeneity at a critical constraint density a cm , where many solution communities form. 
These solution communities serve as precursors for the splitting of the solution space into many solution clusters at a 
larger threshold value ad of constraint density. This work brings a refined picture on how ergodicity of the solution 
space of a CSP finally breaks as the constraint density increases. 

In spin-glass models with multi-spin interactions, the control parameter is often the temperature. The method 
presented here can also be used to study how the configuration spaces of these systems evolve with temperature. We 
suggest that similar heterogeneity transitions will occur before the clustering (or dynamical) transition. The following 
scenario is expected (see Fig. [4]): at high temperatures the configuration space of a spin-glass or a lattice glass model 
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system is in a homogeneous phase; as the temperature T decreases to certain critical value T cm , many communities of 
configurations form in the configuration space, and the configuration space is then in a heterogeneous but still ergodic 
phase; as T decreases further to T d , the different configuration communities separate into different Gibbs states, and 
the configuration space is no longer ergodic. The values of T cm for the random iGSAT problem and the random 
iGXORSAT problem as a function of the constraint density a will be calculated in a forthcoming publication. A 
related study was reported by Krzakala and Zdeb orova recently on the the adiabatic evolution of single Gibbs states 
of a spin-glass system as a function of temperatur e! 24 * 25 !. 

As the solution space of a CSP or the configuration space of a spin-glass or lattice glass system becomes hetero- 
geneous and the configurations aggregate into many different communities, a stochastic search process based only 
on local rules (e.g., solution space random walking^) or a local dynamical process (e.g., single-particle heat-bath 
dynamics of a lattice glass^) may get slowing down considerably and show heterogeneous behavior. The configuration 
space heterogeneity discussed in thi s pap er probably is deeply connected to the phenomenon of spatial dynamical 
heterogeneity of glass-forming liquida^^D. This research direction will be pursued in future work. 
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